#ifndef MESH_H
#define MESH_H

#include "Domain.h"
#include "Dofs.h"

/**
 * @brief 有限元计算网格类，所有网格从这个类派生；
 *        功能：设置和返回网格大小，计算并返回指定全局编号对应的自由度信息，
 *             返回指定单元的所有自由度的全局编号，以及其他网格信息；
 */
template<unsigned int DIM>
class Mesh
{
public:
    /** 
     * 默认构造函数
     * 
     */
    Mesh() = default;
    /**
     * @brief 该网格中单元个数；
     * 
     * @return long int 
     */
    virtual long int n_element() = 0;
    /**
     * @brief 该网格中所有自由度总数;
     * 
     * @return long int 
     */
    virtual long int n_dofs() = 0;
    /** 
     * @brief x或y坐标轴上划分的段数
     * 
     * @param i i=0返回x上划分的段数，i=1返回y上划分的段数。
     * 
     * @return 
     */
    virtual long int& get_seg(int i) = 0;
    /**
     * @brief x坐标轴上每一段长度；
     * 
     * @return long& 
     */
    virtual double x_h() = 0;
    /**
     * @brief y坐标轴上每一段长度；
     * 
     * @return long& 
     */
    virtual double y_h() = 0;
    /**
     * @brief 返回全局编号为i的自由度
     * 
     * @param i 指定全局编号
     * @return Dofs<DIM> 对应的自由度；
     */
    virtual Dofs<DIM> DofsofIndex(int i) = 0;
    /**
     * @brief 返回指定单元所有全局编号；
     * 
     * @param i 指定单元序数；
     * @return std::vector<int> 存放该单元所有全局编号的数组
     */
    virtual std::vector<int> NodeofEle(int i) = 0;
    /**
     * @brief 返回所有边界上点的全局编号；
     * 
     * @return std::vector<int>& 存放该网格所有边界上自由度全局编号的数组
     */
    virtual std::vector<int>& Boundary() = 0;
    /**
     * @brief 返回计算区域的指针
     * 
     * @return RectangleDomain* 计算区域的指针
     */
    virtual RectangleDomain* RDomain() = 0;
    /**
     * @brief 二维度情况下，第i,j号自由度对应的全局编号
     * 
     * @param i 下标１
     * @param j 下标２
     * @return int 全局编号
     */
    virtual int CortoIdx(int i, int j) = 0;
    /**
     * @brief 判断该index序号点是否一个边界点； 
     * 
     * @param index 需要判断点的全局编号
     * @return true 是边界点
     * @return false 不是边界点
     */
    virtual bool IsBoundary(int index) = 0;
    /** 
     * 析构函数
     * 
     */
    virtual ~Mesh(){delete _RDomain;};

protected:
    RectangleDomain* _RDomain;/*<计算区域的指针*/
    long int _seg[2];/*< _seg[0]和_seg[1]分别表示x和y方向上划分的区间数*/
    std::vector<int> _BndIndex;/*<边界上点的全局编号*/
    /** 
     * 构造_BndIndex;
     * 
     */
    virtual void ReBoundary() = 0;
};

/**
 * @brief P1网格，对应单元为P1单元；
 * 
 */
/**********************************************************/
/*********************** P1 2D ****************************/
/**********************************************************/

/* 如何划分单元格？                     
 * 首先将实际计算区域划分为规则的矩形网格，对于每个矩形网格，
 * 将其第0个和第2个顶点连线，划分为两个三角形，从左至右排序。
 * 2  ——  3 
 * |     /|    如图左所示，0-3-2, 0-1-3 构成两个 P1 单元片，
 * |    / |    且 0-3-2 编号在前。
 * |   /  |
 * |  /   |
 * | /    |
 * |/     |
 * 0  ——  1
 *
 */

class P1Mesh:public Mesh<2>
{

public:
    P1Mesh() = default;
    /** 
     * 根据给定区域和每一维要划分的区间数构建网格；
     * 
     * @param Domain 区域
     * @param seg x,y方向上划分的区间数
     * 
     */
    P1Mesh(RectangleDomain* Domain, std::vector<long int> seg);
    ~P1Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_dofs();
    virtual long int& get_seg(int i);
    virtual double x_h();
    virtual double y_h();
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(_seg[0]+1) + j);};
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
     virtual void ReBoundary();
};



P1Mesh::P1Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();

}



long int P1Mesh::n_element()
{
    return static_cast<long int>(2 * _seg[0] * _seg[1]);
};

long int P1Mesh::n_dofs()
{
    return static_cast<long int>((_seg[0] + 1) * (_seg[1] + 1)) ;
};

long int &P1Mesh::get_seg(int i)
{
    return _seg[i];
};

double P1Mesh::x_h()
{
  return fabs((*_RDomain)[0][0] - (*_RDomain)[1][0]) /
         static_cast<double>(_seg[0]);
};

double P1Mesh::y_h()
{
  return fabs((*_RDomain)[0][1] - (*_RDomain)[3][1]) /
         static_cast<double>(_seg[1]);
};

Dofs<2> P1Mesh::DofsofIndex(int index)
{
    Dofs<2> temp_dofs;
    int i = index % (_seg[0] + 1);
    int j = index / static_cast<int>(_seg[0] + 1);
    double x_c = (*_RDomain)[0][0] + i * x_h();
    double y_c = (*_RDomain)[0][1] + j * y_h();
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    return temp_dofs;
}


std::vector<int> P1Mesh::NodeofEle(int i)
{
    std::vector<int> Idx;
     int a = i / (2 * _seg[0]);
    int b = i % (2 * _seg[1]);
    if (b % 2 == 0)
    {
	Idx.push_back(a * (_seg[0] + 1) + b / 2);
	Idx.push_back(Idx[0] + _seg[0] + 2);
	Idx.push_back(Idx[1] - 1);
    }
    else
    {
	Idx.push_back(a * (_seg[0] + 1) + (b -1) / 2);
	Idx.push_back(Idx[0] + 1);
	Idx.push_back(Idx[1] + _seg[0] + 1);
    }
    return Idx;
}

void P1Mesh::ReBoundary()
{
    for(int i = 0;i <= _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_dofs()- i - 1);
    }
    for(int j = 1;j <= _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j*(_seg[0] + 1));
        _BndIndex.push_back((j + 1)*(_seg[0] + 1) - 1);
    }
}
/**
 * @brief P2网格，对应单元为P2单元；
 * 
 */
/**********************************************************/
/*********************** P2 2D ****************************/
/**********************************************************/

/* 如何划分单元格？                     
 * 首先将实际计算区域划分为规则的矩形网格，对于每个矩形网格，
 * 将其第0个和第8个顶点连线，划分为两个三角形，从左至右排序。
 * 从左下角为第一个号，逆时针排序．
 * 6  —7—  8 
 * |     /|    如图左所示，0-8-6-4-7-3, 0-2-8-1-5-4 构成两个 P2 单元片，
 * |    / |    且 0-8-6-4-7-3 编号在前。
 * 3   4  5
 * |  /   |
 * | /    |
 * |/     |
 * 0  —1—  2
 *
 */

class P2Mesh:public Mesh<2>
{
public:
    P2Mesh(){};
    /** 
     * 根据给定区域和每一维要划分的区间数构建网格；
     * 
     * @param Domain 区域
     * @param seg x,y方向上划分的区间数
     * 
     */
    P2Mesh(RectangleDomain* Domain, std::vector<long int> seg);
    ~P2Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_dofs();
    virtual long int& get_seg(int i){return _seg[i];};
    virtual double x_h() {
	return fabs((*_RDomain)[0][0] - (*_RDomain)[1][0]) / static_cast<double>(_seg[0]);
    };
    virtual double y_h() {
	return fabs((*_RDomain)[0][1] - (*_RDomain)[3][1]) / static_cast<double>(_seg[1]);      
    };
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(2*_seg[0]+1) + j);};
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
     virtual void ReBoundary();
};



P2Mesh::P2Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();

}



long int P2Mesh::n_element()
{
    return static_cast<long int>(2 * _seg[0] * _seg[1]);
}

long int P2Mesh::n_dofs()
{
    return static_cast<long int>((2 * _seg[0] + 1) * (2 * _seg[1] + 1)) ;
}


Dofs<2> P2Mesh::DofsofIndex(int index)
{
    Dofs<2> temp_dofs;
    int i = index % (2 * _seg[0] + 1);
    int j = index / static_cast<int>(2 * _seg[0] + 1);
    double x_c = (*_RDomain)[0][0] + i * x_h() / 2.0;
    double y_c = (*_RDomain)[0][1] + j * y_h() / 2.0;
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    return temp_dofs;
}


std::vector<int> P2Mesh::NodeofEle(int i)
{
    std::vector<int> Idx;
    int a = i / (2 * _seg[0]);
    int b = i % (2 * _seg[0]);
    if (b % 2 == 0)
    {
	Idx.push_back(a * 2 * (2 * _seg[0] + 1) + b);
	Idx.push_back(Idx[0] + 2 * (2 * _seg[0] + 1) + 2);
	Idx.push_back(Idx[1] - 2);
        Idx.push_back(Idx[0] + (2 * _seg[0] + 1) + 1);
        Idx.push_back(Idx[1] - 1);
        Idx.push_back(Idx[3] - 1);
    }
    else
    {
	Idx.push_back(a * 2 *(2 * _seg[0] + 1) + b - 1);
	Idx.push_back(Idx[0] + 2);
	Idx.push_back(Idx[1] + 2 * (2 * _seg[0] + 1));
        Idx.push_back(Idx[0] + 1);
        Idx.push_back(Idx[1] + 2 * _seg[0] + 1);
        Idx.push_back(Idx[4] - 1);
    }
    return Idx;
}

void P2Mesh::ReBoundary()
{
    for(int i = 0;i <= 2 * _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_dofs()- i - 1);
    }
    for(int j = 1;j <= 2 * _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j* (2  * _seg[0] + 1));
        _BndIndex.push_back((j + 1)*(2 * _seg[0] + 1) - 1);
    }
}


/**********************************************************/
/*********************** P3 2D ****************************/
/**********************************************************/
/**
 * @brief P3网格，对应单元为P3单元；
 * 
 */
class P3Mesh:public Mesh<2>
{
public:
    P3Mesh(){};
    /** 
     * 根据给定区域和每一维要划分的区间数构建网格；
     * 
     * @param Domain 区域
     * @param seg x,y方向上划分的区间数
     * 
     */
    P3Mesh(RectangleDomain* Domain, std::vector<long int> seg);
    ~P3Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_dofs();
    virtual long int& get_seg(int i){return _seg[i];};
    virtual double x_h() {
	return fabs((*_RDomain)[0][0] - (*_RDomain)[1][0]) / static_cast<double>(_seg[0]);
    };
    virtual double y_h() {
	return fabs((*_RDomain)[0][1] - (*_RDomain)[3][1]) / static_cast<double>(_seg[1]);      
    };
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(3*_seg[0]+1) + j);};
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
     virtual void ReBoundary();
};



P3Mesh::P3Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();

}



long int P3Mesh::n_element()
{
    return static_cast<long int>(2 * _seg[0] * _seg[1]);
}

long int P3Mesh::n_dofs()
{
    return static_cast<long int>((3 * _seg[0] + 1) * (3 * _seg[1] + 1)) ;
}


Dofs<2> P3Mesh::DofsofIndex(int index)
{
    Dofs<2> temp_dofs;
    int i = index % (3 * _seg[0] + 1);
    int j = index / static_cast<int>(3 * _seg[0] + 1);
    double x_c = (*_RDomain)[0][0] + i * x_h() / 3.0;
    double y_c = (*_RDomain)[0][1] + j * y_h() / 3.0;
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    return temp_dofs;
}

std::vector<int> P3Mesh::NodeofEle(int i)
{
    std::vector<int> Idx;
    int a = i / (2 * _seg[0]);
    int b = i % (2 * _seg[0]);
    if (b % 2 == 0)
    {
	Idx.push_back(a * 3 * (3 * _seg[0] + 1) + 3*b/2);
	Idx.push_back(Idx[0] + 3 * (3 * _seg[0] + 1) + 3);
	Idx.push_back(Idx[1] - 3);
	Idx.push_back(Idx[0] + (3 * _seg[0] + 1) + 1);
	Idx.push_back(Idx[0] + 2 * (3 * _seg[0] + 1) + 2);	
	Idx.push_back(Idx[1] - 1);
	Idx.push_back(Idx[1] - 2);
	Idx.push_back(Idx[0] + 2 * (3 * _seg[0] + 1));
	Idx.push_back(Idx[0] + (3 * _seg[0] + 1));
	Idx.push_back(Idx[4] - 1);
    }
    else
    {
	Idx.push_back(a * 3 *(3 * _seg[0] + 1) + 3*(b-1)/2);
	Idx.push_back(Idx[0] + 3);
	Idx.push_back(Idx[1] + 3 * (3 * _seg[0] + 1));
        Idx.push_back(Idx[0] + 1);
	Idx.push_back(Idx[0] + 2);
        Idx.push_back(Idx[1] + (3 * _seg[0] + 1));
	Idx.push_back(Idx[1] + 2 * (3 * _seg[0] + 1));
        Idx.push_back(Idx[6] - 1);
	Idx.push_back(Idx[5] - 2);
	Idx.push_back(Idx[5] - 1);
    }
    return Idx;
}

void P3Mesh::ReBoundary()
{
    for(int i = 0;i <= 3 * _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_dofs()- i - 1);
    }
    for(int j = 1;j <= 3 * _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j* (3  * _seg[0] + 1));
        _BndIndex.push_back((j + 1)*(3 * _seg[0] + 1) - 1);
    }
}

/**********************************************************/
/*********************** P4 2D ****************************/
/**********************************************************/
/**
 * @brief P4网格，对应单元为P4单元；
 * 
 */
class P4Mesh:public Mesh<2>
{
public:
    P4Mesh(){};
    /** 
     * 根据给定区域和每一维要划分的区间数构建网格；
     * 
     * @param Domain 区域
     * @param seg x,y方向上划分的区间数
     * 
     */
    P4Mesh(RectangleDomain* Domain, std::vector<long int> seg);
    ~P4Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_dofs();
    virtual long int& get_seg(int i){return _seg[i];};
    virtual double x_h() {
	return fabs((*_RDomain)[0][0] - (*_RDomain)[1][0]) / static_cast<double>(_seg[0]);
    };
    virtual double y_h() {
	return fabs((*_RDomain)[0][1] - (*_RDomain)[3][1]) / static_cast<double>(_seg[1]);      
    };
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(4*_seg[0]+1) + j);};
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
     virtual void ReBoundary();
};



P4Mesh::P4Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();

}



long int P4Mesh::n_element()
{
    return static_cast<long int>(2 * _seg[0] * _seg[1]);
}

long int P4Mesh::n_dofs()
{
    return static_cast<long int>((4 * _seg[0] + 1) * (4 * _seg[1] + 1)) ;
}


Dofs<2> P4Mesh::DofsofIndex(int index)
{
    Dofs<2> temp_dofs;
    int i = index % (4 * _seg[0] + 1);
    int j = index / static_cast<int>(4 * _seg[0] + 1);
    double x_c = (*_RDomain)[0][0] + i * x_h() / 4.0;
    double y_c = (*_RDomain)[0][1] + j * y_h() / 4.0;
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    return temp_dofs;
}

std::vector<int> P4Mesh::NodeofEle(int i)
{
    std::vector<int> Idx;
    int a = i / (2 * _seg[0]);
    int b = i % (2 * _seg[0]);
    if (b % 2 == 0)
    {
	Idx.push_back(a * 4 * (4 * _seg[0] + 1) + 2*b);
	Idx.push_back(Idx[0] + 4 * (4 * _seg[0] + 1) + 4);
	Idx.push_back(Idx[1] - 4);
	Idx.push_back(Idx[0] + (4 * _seg[0] + 1) + 1);
	Idx.push_back(Idx[0] + 2 * (4 * _seg[0] + 1) + 2);
	Idx.push_back(Idx[0] + 3 * (4 * _seg[0] + 1) + 3);
	Idx.push_back(Idx[1] - 1);
	Idx.push_back(Idx[1] - 2);
	Idx.push_back(Idx[1] - 3);
	Idx.push_back(Idx[0] + 3 * (4 * _seg[0] + 1));
	Idx.push_back(Idx[0] + 2 * (4 * _seg[0] + 1));
	Idx.push_back(Idx[0] + (4 * _seg[0] + 1));
	Idx.push_back(Idx[4] - 1);
	Idx.push_back(Idx[5] - 1);
	Idx.push_back(Idx[5] - 2);
    }
    else
    {
	Idx.push_back(a * 4 *(4 * _seg[0] + 1) + 2*(b-1));
	Idx.push_back(Idx[0] + 4);
	Idx.push_back(Idx[1] + 4 * (4 * _seg[0] + 1));
        Idx.push_back(Idx[0] + 1);
	Idx.push_back(Idx[0] + 2);
	Idx.push_back(Idx[0] + 3);
        Idx.push_back(Idx[1] + (4 * _seg[0] + 1));
	Idx.push_back(Idx[1] + 2 * (4 * _seg[0] + 1));
	Idx.push_back(Idx[1] + 3 * (4 * _seg[0] + 1));
        Idx.push_back(Idx[8] - 1);
	Idx.push_back(Idx[7] - 2);
	Idx.push_back(Idx[6] - 3);
	Idx.push_back(Idx[6] - 2);
	Idx.push_back(Idx[6] - 1);
	Idx.push_back(Idx[7] - 1);
    }
    return Idx;
}

void P4Mesh::ReBoundary()
{
    for(int i = 0;i <= 4 * _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_dofs()- i - 1);
    }
    for(int j = 1;j <= 4 * _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j* (4  * _seg[0] + 1));
        _BndIndex.push_back((j + 1)*(4 * _seg[0] + 1) - 1);
    }
}



/**
 * @brief Q1网格，对应单元为Q1单元；
 * 
 */
/**********************************************************/
/*********************** Q1 2D ****************************/
/**********************************************************/

class Q1Mesh:public Mesh<2>
{
public:
    Q1Mesh(){};
    Q1Mesh(RectangleDomain*,std::vector<long int>);
    ~Q1Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_dofs();
    virtual long int& get_seg(int i); 
    virtual double x_h();
    virtual double y_h();
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(_seg[0]+1) + j);};
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
    virtual void ReBoundary();
};




Q1Mesh::Q1Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif

}


long int Q1Mesh::n_element()
{
    return static_cast<long int>(_seg[0] * _seg[1]);
}

long int Q1Mesh::n_dofs()
{
    return static_cast<long int>((_seg[0] + 1) * (_seg[1] + 1));
}
long int& Q1Mesh::get_seg(int i)
{
    return _seg[i];
};
double Q1Mesh::x_h()
{
  return fabs((*_RDomain)[0][0] - (*_RDomain)[1][0]) /
         static_cast<double>(_seg[0]);
};

double Q1Mesh::y_h()
{
  return fabs((*_RDomain)[0][1] - (*_RDomain)[3][1]) /
         static_cast<double>(_seg[1]);
};



/** 
 * @brief 返回全局编号对应的自由度
 * 
 * @param index 
 * 
 * @return 
 */
Dofs<2> Q1Mesh::DofsofIndex(int index)
{
    Dofs<2> temp_dofs;
    int i = index % (_seg[0] + 1);
    int j = index / static_cast<int>(_seg[0] + 1);
    double x_c = (*_RDomain)[0][0] + i * x_h();
    double y_c = (*_RDomain)[0][1] + j * y_h();
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    return temp_dofs;
} 
 
/** 
 * @brief 按顺序返回指定单元包含的所有自由度的编号
 * 
 * @param i 
 * 
 * @return 
 */
std::vector<int> Q1Mesh::NodeofEle(int i)
{
    std::vector<int> Idx;
    int a = i / (_seg[0]);
    int b = i % (_seg[0]);
    Idx.push_back(a * (_seg[0] + 1) + b);
    Idx.push_back(Idx[0] + 1);
    Idx.push_back(Idx[1] + (_seg[0] + 1));
    Idx.push_back(Idx[2] - 1);
    return Idx;
}
void Q1Mesh::ReBoundary()
{
    for(int i = 0;i <= _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_dofs()- i - 1);
    }
    for(int j = 1;j <= _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j*(_seg[0] + 1));
        _BndIndex.push_back((j + 1)*(_seg[0] + 1) - 1);
    }
}



/**********************************************************/
/*********************** Q2 2D ****************************/
/**********************************************************/
class Q2Mesh:public Mesh<2>
{
public:

    Q2Mesh(){};
    Q2Mesh(RectangleDomain*,std::vector<long int>);
    ~Q2Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_dofs();
    virtual long int& get_seg(int i);
    virtual double x_h();
    virtual double y_h();
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(2*_seg[0]+1) + j);};//Q2 每行有两倍多的点。
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
    virtual void ReBoundary();
private:
    RectangleDomain* _RDomain;
    long int _seg[2];
    std::vector<int> _BndIndex;

};




Q2Mesh::Q2Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif

}

long int Q2Mesh::n_element()
{
    return static_cast<long int>(_seg[0] * _seg[1]);
}

long int Q2Mesh::n_dofs()
{
    return static_cast<long int>((2 * _seg[0] + 1) * (2 * _seg[1] + 1));
}

long int& Q2Mesh::get_seg(int i)
{
    return _seg[i];
};
double Q2Mesh::x_h()
{
  return fabs((*_RDomain)[0][0] - (*_RDomain)[1][0]) /
         static_cast<double>(_seg[0]);
};

double Q2Mesh::y_h()
{
  return fabs((*_RDomain)[0][1] - (*_RDomain)[3][1]) /
         static_cast<double>(_seg[1]);
};



/** 
 * @brief 返回全局编号对应的自由度
 * 
 * @param index 
 * 
 * @return 
 */
Dofs<2> Q2Mesh::DofsofIndex(int index)
{

    Dofs<2> temp_dofs;
    int i = index % (2 * _seg[0] + 1);
    int j = index / static_cast<int>(2 * _seg[0] + 1);
    double x_c = (*_RDomain)[0][0] + i * x_h() /2.0;
    double y_c = (*_RDomain)[0][1] + j * y_h() /2.0;
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    return temp_dofs;
} 

/** 
 * @brief 按顺序返回指定单元包含的所有自由度的编号
 * 
 * @param i 
 * 
 * @return 
 */
std::vector<int> Q2Mesh::NodeofEle(int i)
{
    std::vector<int> Idx;
    int a = i / (_seg[0]);
    int b = i % (_seg[0]);
    Idx.push_back(a * 2 * (2 * _seg[0] + 1) + 2 * b);
    Idx.push_back(Idx[0] + 2);
    Idx.push_back(Idx[1] + 2 * (2 * _seg[0] + 1));
    Idx.push_back(Idx[2] - 2);
    Idx.push_back(Idx[0] + 1);
    Idx.push_back(Idx[1] + (2 * _seg[0] + 1));
    Idx.push_back(Idx[2] - 1);
    Idx.push_back(Idx[0] + (2 * _seg[0] + 1));
    Idx.push_back(Idx[7] + 1);
    return Idx;
}
void Q2Mesh::ReBoundary()
{
    for(int i = 0;i <= 2 * _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_dofs()- i - 1);
    }
    for(int j = 1;j <= 2 * _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j* (2  * _seg[0] + 1));
        _BndIndex.push_back((j + 1)*(2 * _seg[0] + 1) - 1);
    }
}

/**
 * @brief Q3网格，对应单元为Q3单元；
 * 
 */
/**********************************************************/
/*********************** Q3 2D ****************************/
/**********************************************************/
class Q3Mesh:public Mesh<2>
{
public:

    Q3Mesh(){};
    Q3Mesh(RectangleDomain*,std::vector<long int>);
    ~Q3Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_dofs();
    virtual long int& get_seg(int i);
    virtual double x_h();
    virtual double y_h();
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(3*_seg[0]+1) + j);};//Q3 每行有三倍多的点。
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
    virtual void ReBoundary();
private:
    RectangleDomain* _RDomain;
    long int _seg[2];
    std::vector<int> _BndIndex;

};




Q3Mesh::Q3Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif

}

long int Q3Mesh::n_element()
{
    return static_cast<long int>(_seg[0] * _seg[1]);
}

long int Q3Mesh::n_dofs()
{
    return static_cast<long int>((3 * _seg[0] + 1) * (3 * _seg[1] + 1));
}

long int& Q3Mesh::get_seg(int i)
{
    return _seg[i];
};
double Q3Mesh::x_h()
{
  return fabs((*_RDomain)[0][0] - (*_RDomain)[1][0]) /
         static_cast<double>(_seg[0]);
};

double Q3Mesh::y_h()
{
  return fabs((*_RDomain)[0][1] - (*_RDomain)[3][1]) /
         static_cast<double>(_seg[1]);
};



/** 
 * @brief 返回全局编号对应的自由度
 * 
 * @param index 
 * 
 * @return 
 */
Dofs<2> Q3Mesh::DofsofIndex(int index)
{

    Dofs<2> temp_dofs;
    int i = index % (3 * _seg[0] + 1);
    int j = index / static_cast<int>(3 * _seg[0] + 1);
    double x_c = (*_RDomain)[0][0] + i * x_h() /3.0;
    double y_c = (*_RDomain)[0][1] + j * y_h() /3.0;
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    return temp_dofs;
} 

/** 
 * @brief 按顺序返回指定单元包含的所有自由度的编号
 * 
 * @param i 
 * 
 * @return 
 */
std::vector<int> Q3Mesh::NodeofEle(int i)
{
    std::vector<int> Idx;
    int a = i / (_seg[0]);
    int b = i % (_seg[0]);
    Idx.push_back(a * 3 * (3 * _seg[0] + 1) + 3 * b);
    Idx.push_back(Idx[0] + 3);
    Idx.push_back(Idx[1] + 3 * (3 * _seg[0] + 1));
    Idx.push_back(Idx[2] - 3);
    Idx.push_back(Idx[0] + 1);
    Idx.push_back(Idx[0] + 2);
    Idx.push_back(Idx[1] + (3 * _seg[0] + 1));
    Idx.push_back(Idx[1] + 2 * (3 * _seg[0] + 1));
    Idx.push_back(Idx[2] - 1);
    Idx.push_back(Idx[2] - 2);
    Idx.push_back(Idx[0] + 2 * (3 * _seg[0] + 1));
    Idx.push_back(Idx[0] + (3 * _seg[0] + 1));
    Idx.push_back(Idx[11] + 1);
    Idx.push_back(Idx[6] - 1);
    Idx.push_back(Idx[7] - 1);
    Idx.push_back(Idx[10] + 1);
    return Idx;
}
void Q3Mesh::ReBoundary()
{
    for(int i = 0;i <= 3 * _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_dofs()- i - 1);
    }
    for(int j = 1;j <= 3 * _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j* (3  * _seg[0] + 1));
        _BndIndex.push_back((j + 1)*(3 * _seg[0] + 1) - 1);
    }
}





/**
 * @brief Q4网格，对应单元为Q4单元；
 * 
 */
/**********************************************************/
/*********************** Q4 2D ****************************/
/**********************************************************/
class Q4Mesh:public Mesh<2>
{
public:

    Q4Mesh(){};
    Q4Mesh(RectangleDomain*,std::vector<long int>);
    ~Q4Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_dofs();
    virtual long int& get_seg(int i);
    virtual double x_h();
    virtual double y_h();
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(4*_seg[0]+1) + j);};//Q4 每行有四倍多的点。
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
    virtual void ReBoundary();
private:
    RectangleDomain* _RDomain;
    long int _seg[2];
    std::vector<int> _BndIndex;

};




Q4Mesh::Q4Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif

}

long int Q4Mesh::n_element()
{
    return static_cast<long int>(_seg[0] * _seg[1]);
}

long int Q4Mesh::n_dofs()
{
    return static_cast<long int>((4 * _seg[0] + 1) * (4 * _seg[1] + 1));
}

long int& Q4Mesh::get_seg(int i)
{
    return _seg[i];
};
double Q4Mesh::x_h()
{
  return fabs((*_RDomain)[0][0] - (*_RDomain)[1][0]) /
         static_cast<double>(_seg[0]);
};

double Q4Mesh::y_h()
{
  return fabs((*_RDomain)[0][1] - (*_RDomain)[3][1]) /
         static_cast<double>(_seg[1]);
};



/** 
 * @brief 返回全局编号对应的自由度
 * 
 * @param index 
 * 
 * @return 
 */
Dofs<2> Q4Mesh::DofsofIndex(int index)
{

    Dofs<2> temp_dofs;
    int i = index % (4 * _seg[0] + 1);
    int j = index / static_cast<int>(4 * _seg[0] + 1);
    double x_c = (*_RDomain)[0][0] + i * x_h() /4.0;
    double y_c = (*_RDomain)[0][1] + j * y_h() /4.0;
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    return temp_dofs;
} 

/** 
 * @brief 按顺序返回指定单元包含的所有自由度的编号
 * 
 * @param i 
 * 
 * @return 
 */
std::vector<int> Q4Mesh::NodeofEle(int i)
{
    std::vector<int> Idx;
    int a = i / (_seg[0]);
    int b = i % (_seg[0]);
    Idx.push_back(a * 4 * (4 * _seg[0] + 1) + 4 * b);
    Idx.push_back(Idx[0] + 4);
    Idx.push_back(Idx[1] + 4 * (4 * _seg[0] + 1));
    Idx.push_back(Idx[2] - 4);
    Idx.push_back(Idx[0] + 1);
    Idx.push_back(Idx[0] + 2);
    Idx.push_back(Idx[0] + 3);
    Idx.push_back(Idx[1] + (4 * _seg[0] + 1));
    Idx.push_back(Idx[1] + 2 * (4 * _seg[0] + 1));
    Idx.push_back(Idx[1] + 3 * (4 * _seg[0] + 1));
    Idx.push_back(Idx[2] - 1);
    Idx.push_back(Idx[2] - 2);
    Idx.push_back(Idx[2] - 3);
    Idx.push_back(Idx[0] + 3 * (4 * _seg[0] + 1));
    Idx.push_back(Idx[0] + 2 * (4 * _seg[0] + 1));
    Idx.push_back(Idx[0] + (4 * _seg[0] + 1));
    Idx.push_back(Idx[7] - 3);
    Idx.push_back(Idx[7] - 2);
    Idx.push_back(Idx[7] - 1);
    Idx.push_back(Idx[8] - 1);
    Idx.push_back(Idx[9] - 1);
    Idx.push_back(Idx[9] - 2);
    Idx.push_back(Idx[9] - 3);
    Idx.push_back(Idx[8] - 3);
    Idx.push_back(Idx[8] - 2);
    return Idx;
}
void Q4Mesh::ReBoundary()
{
    for(int i = 0;i <= 4 * _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_dofs()- i - 1);
    }
    for(int j = 1;j <= 4 * _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j* (4  * _seg[0] + 1));
        _BndIndex.push_back((j + 1)*(4 * _seg[0] + 1) - 1);
    }
}


#endif /* MESH_H */
